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For large but finite systems the static properties of the infinite ranged Sherrington-Kirkpatrick 
model are numerically investigated in the entire the glass regime. The approach is based on the 
modified Thouless-Anderson-Palmer equations in combination with a phenomenological relaxational 
dynamics used as a numerical tool. For all temperatures and all bond configurations stable and meta 
stable states are found. Following a discussion of the finite size effects, the static properties of the 
state of lowest free energy are presented in the presence of a homogeneous magnetic field for all 
temperatures below the spin glass temperature. Moreover some characteristic features of the meta 
stable states are presented. These states exist in finite temperature intervals and disappear via local 
saddle node bifurcations. Numerical evidence is found that the excess free energy of the meta stable 
states remains finite in the thermodynamic limit. This implies a the ' mult i- valley' structure of the 
free energy on a sub-extensive scale. 

PACS numbers: 75.10.Nr, 05.50.+q, 87.10.+C 



I. INTRODUCTION 

The Sherrington-Kirkpatrick (SK) model Q for Ising 
spins with quenched random bonds is the simplest rep- 
resentative of a class of long-ranged models all describ- 
ing successfully the interesting phenomena of spin glasses 
^ [|. In addition to this success in physical questions, 
the research on these models has been fruitful and stim- 
ulating in optimization problems, in the understanding 
of neural networks and in the field of learning machines 

Two complimentary but conceptually different ap- 
proaches exist for the analysis of these models. The first 
approach is the replica method applied already by Sher- 
rington and Kirkpatrick to their model in the original 
work. After breaking the replica symmetry it was Parisi 
[^| who first found a satisfactory solution of the SK model 
in the spin glass regime. 

Within the replica theory and the related dynamical 
approaches Q , an averaging over the quenched variables 
is performed and therefore only the calculation of bond- 
averaged quantities is possible. All problems that de- 
pend on a special, but arbitrary configuration of the ran- 
dom bonds cannot be solved within the framework of the 
replica theory. Such problems are the generic case in the 
nonphysical applications named above. In optimization 
problems, for instance, the main question is to find the 
optimal solution for a fixed configuration. 

The complimentary approach of Thouless-Anderson- 
Palmer (TAP) to investigate spin glass models does 
not perform the bond average and permits the treatment 
of problems depending on specific configurations. For 
other questions which are expected to be independent of 
the special configuration, such as all macroscopic phys- 
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ical quantities, self averaging happens. This is due to 
the fact that the random interaction matrices have well 
known asymptotic properties in the thermodynamic limit 
S. The situation is in principle similar to the central 
limit theorem in probability theory, where large numbers 
of random variables also permit the calculation of macro- 
scopic quantities which hold for nearly every realization 
of the random variables. Thus the investigation of one 
or some representative systems is sufficient and the bond 
average is not needed. 

The TAP equations have been well established for more 
than two decades and several alternative derivations are 
known ]|, ||, ^| . Nevertheless the TAP approach is still a 
field of current interest. This is due to the importance of 
the approach to numerous interesting problems. More- 
over it is suspected that not all aspects of this approach 
have yet been worked out. Very recently the author |T(]] 
has reanalyzed the stability of the TAP equations. He 
concluded that the reason for the breakdown at the spin 
glass instability is an inconsistency in the value of the lo- 
cal susceptibility and showed that unstable states cannot 
be described by the original TAP equations. Remov- 
ing these inconsistencies by a self-consistent treatment, 
modified TAP equations below the instability have been 
presented [l0| . 

One of the pressing questions is the characteristic fea- 
tures of the pure states of the SK model, which are the 
stable solutions of the TAP equations. It is obvious that 
these states are essential for a direct understanding of the 
SK model. Moreover the pure states are of importance 
for the interpretation of the more formal results of the 
replica theory |2| . 

Apart from the research JTT[ , showing that the number 
of the pure states is very large at low temperatures, two 
approaches exist to investigate the characteristic features 
of the pure states. Both the work of Bray and Moore Jl^] 
and the work of Nemoto and Takayama 1 1 3 try to find 



explicitly numerical solutions of the TAP equations for 
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large but finite systems. Only in rare cases (about 15 
% of the systems with different bond configurations) do 
these investigations yield solutions of the TAP equations 
which satisfy the validity condition. To treat the other 
cases Nemoto and Takayama have developed an interest- 
ing but indirect method to obtain at least approximative 
solutions which are expected to converge to true solu- 
tions in the thermodynamic limit. The applied numeri- 
cal method of |l3j did not work well for a cooling of the 
sample and therefore results for this case have not been 
published. 

In this work we will demonstrate that all the known 
difficulties to find numerically the pure states disappear 
if the modified TAP equations (Hj are applied instead 
of the original ones. For all configurations of the bonds, 
stable solutions of these TAP equations will be found for 
finite systems employing a simple relaxational dynamics. 
There is no further need to apply the indirect method of 
p3| and our approach permits the sample to be cooled 
down, as is shown to be essential to find the low lying 
states. In addition to these general goals some charac- 
teristic features of pure states will be worked out leading 
to a better understanding of the SK spin glass. 

In Section II the basics of the modified TAP equations, 
some essentials of the applied dynamics and the numeri- 
cal procedure are presented. A main subject of section III 
is the investigation of finite size effects which are found 
to be of some importance. The numerical results of the 
state of lowest free energy for the relevant static quanti- 
ties will be given for all temperatures below the critical 
spin glass temperature. Next we focus on the meta stable 
states and discuss some characteristic features. Finally, 
the concluding remarks can be found in section IV. 



II. BASIC EQUATIONS AND METHOD 

A. The modified TAP equations 

The Hamiltonian of the SK model of a system of N 
Ising spins [Si = ±1) in the presence of local external 
fields hi is described by 



H — y'^^JijSiSj ^ hjSj 



(1) 



The bonds are independent random variables with 
zero means and standard deviations N~ x / 2 (where the 
latter scaling fixes the spin glass temperature to T = 1 ) . 

The modified TAP mean field equations (l(J] describe 
the statics of the system (||) in the thermodynamic limit. 
They are given by the set of equations for the local mag- 
netizations , with rrii = (Si)p 



rrii = tanh /3{/ij + Ji- 
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together with the local susceptibility 
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and with 

T = 
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x = 1-/3 2 (1- 2q 2 + q 4 ) 



forz<0 (5) 
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(7) 



was introduced. 

The condition x — re pre sents the central spin glass 
instability condition |§|, |l4f. Above the instability the 

local susceptibility xi reduces to the isothermal value \ 



X+ = xP } = 0(1 - 92) 



for x > 



(8) 



which is in complete agreement with the original TAP ap- 
proach ||, 0, [|. Essential differences, however, result 
below the instability as 



r- > o , 



Xi 7^ Xi for a; < 

holds for the modified TAP equations. According to |ITJ 
the quantity T is proportional to the density of zero eigen- 
values of the inverse susceptibility matrix. Thus all the 
states with x < are unstable. As thermodynamics does 
not apply to unstable states, a conflict with the exact 
thermodynamic relation xi — Xi does not exist. 

As worked out in jlo| the modified Tap equations re- 
sult from the simple fact that the eqs.(||) together with 
the definition xi = N~ 1 J2 i drrii/dhi are complete and 
sufficient to determine the m,. From the strict thermo- 
dynamic point of view the modified equations are equiv- 
alent to the original TAP equations as thermodynamics 
is a priori limited to stable states. Note, however, that 
unstable states are essential for discussions based on free 
energy landscapes or for the extension to dynamics of the 
present approach. For a further discussion of this point 
we refer to [^TJ and to the end of the next subsection. 

For later use we give the well known expressions (|, |j, 
0, ^] for the energy 
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and the entropy 



S= - 
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from which the free energy 

F = U - TS 



(12) 



can be calculated. Note that these expressions were ob- 
tained by thermodynamic approaches. Hence they are 
restricted to stable states. 



B. Glauber dynamics 

As in ||, [l^] , Glauber dynamics in mean field approx- 
imation will phenomenologically be added to the model. 
Measuring the time t in units of the relaxation time, the 
purely relaxational equations of motion are given by 

rhi(t) = -mi + t&nh/3{hi + 2J Jymj - miXi(t)}. 

(13) 

The local susceptibility Xz(i) is related to the rrii(t) via 
eq.(|^) and via eq.(Q) or eq.(gj) depending on the instan- 
taneous value of x(t). The fixpoints of these equations 
of motion obviously coincide with the solutions of the 
modified TAP equations. 

To analyze the dynamical stability of these fixpoints 
the equations of motion are linearized near a fixpoint 
rrii. Setting 



8rrii(t) — rrii(t) — nii 



(14) 



and neglecting as usual the N 1 order terms oc d\i /drrii , 
the linearized equations take in matrix notation the form 



5m = — XoX l 5m 



(15) 



where Xo represents the static susceptibility matrix for 
noninteracting Ising spins 



and where 



(Xo)y = 0(1 - k 



X 1 =X 1 +»1-J 



(16) 



(17) 



is the inverse of the susceptibility matrix of the modi- 
fied TAP approach jn^. For thermodynamic stable or 
marginal stable states x~ l is positive semi-definite. 

According to eq.([l5|) the dynamical stability is gov- 
erned the matrix A = XoX -1 - Let I be an eigenvector of 
A satisfying the eigenvalue equation A I — XL Multipli- 
cation with Xo 1 leads to x~ X 1 — A Xq 1 I anc ^ to 



A 



I 



I 



l-Xa'l 



Immediately A > results from the dcfinitcncss of xo an d 
of x 1 f° r stable (or meta stable) static states. Static 
stability therefore implies dynamical stability. By a sim- 
ilar argument, the statement can be proven in the op- 
posite direction, leading to a one-to-one correspondence 
between the two types of stability. 



The latter result implies that the stable solutions of 
the TAP equations can be found by integration of the 
equations of motions (13). The system relaxes to a sta- 
ble fixpoint which is certainly located outside the regime 
x < 0. For transient times, however, the system can en- 
ter the regime x(t) < and thus the flow of eqs. (fl3|) is 
needed in both regions x(t) < and x(t) > 0. These ar- 
guments demonstrate the relevance of the unstable states 
and consequently the need for the modified TAP equa- 
tions for the present approach. It should be added that 
equations of motion of the above type but based on the 
original TAP equations usually lead to incorrect results 
in the spin glass region. For this case the system gener- 
ally relaxes to paramagnetic solutions. 



C. Systems of finite size 

The modified TAP equations are exact in the thermo- 
dynamic limit but are approximative for systems of finite 
size. As numerical studies can only be performed for fi- 
nite systems the two approximations yielding to eqs.(^|j^) 
for finite N are recalled. 

First of all eqs.(^) holds only to leading order in N^ 1 
as sub-extensive contributions have been neglected (com- 
pare e.g. appendix of |l(J). Further corrections to eq (|f ) 
result from the application jnj of the Pastur theorem |lq] 
which is again only exact in the thermodynamic limit. 
The explicit form of these two types of corrections are 
not known and seem to be hard to determine analyti- 
cally. Hence the finite size effects resulting from the two 
approximations will be investigated empirically in the nu- 
merical section. 

Due to the sub-extensive corrections the border of sta- 
bility x — is expected to be only approximative for 
finite systems. This implies that unstable solutions with 
x > and stable solutions with x < may exist for sys- 
tems of finite size. Such stable solutions with x < and 
T > are indeed found in the following. 

Finally it is pointed out that the results of subsection 
B are not restricted to the thermodynamic limit and hold 
for systems of all sizes. Thus the fixpoints obtained by 
integration of the eqs.([l3|) always correspond to stable 
solutions of the modified TAP equations. 



D. Numerical procedure 

The numerical investigations are performed for systems 
up to a size of N = 225 with binary distributions of the 
bonds Jij — ±A -1 / 2 using the standard routine 'Nlnte- 
grate' of Mathematica on workstations. 

The integration of the equations of motion ([l3]) in 
the regime x(t) > is straightforward. In the region 
x(t) < 0, however, a special procedure is employed and 
(/3r) 2 is treated as additional dynamical variable. The 
time derivative of eq.([5]) yields the necessary additional 
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FIG. 1: Approach to the spin glass instability from 
below: T 2 versus 7V~ 2 / 3 at the temperature T=.5. The data 
points represent the averages over about 85, over 28 and over 
25 samples for N < 121, for N = 196 and for N = 225, 
respectively. The length of each error bar is two standard 
deviations. 



FIG. 2: Approach to the spin glass instability from 
below: x — 1 — /3 2 (1 — 2q? +94) versus N~ 2 ^ 3 at the tempera- 
ture T=.5. The data points represent the averages over about 
85, over 28 and over 25 samples for N < 121, for N = 196 
and for N = 225, respectively. The length of each error bar 
is two standard deviations. 



equation of motion and the initial value of (/JT) 2 is de- 
termined from eqs.(^|) for t — which finally guaranties 
that eq.(|5|) is satisfied for all times. 

During the dynamic evolution x(t) generally changes 
sign, which implies a change of the numerical treat- 
ment. Therefore an alternative, approximative approach 
is partly applied. Instead of eq.(||) the equation 



N ^ 1 



/3 2 (l-m 2 ) 2 



NT 2 ' l + r 2 /3 2 (l-m 2 ) 2 



(18) 



is used, where c is a small positive constant with typical 
values from 10~ 3 to 10 -5 . With the supplementary term 
c/(NT 2 ) eq.(Q has a solution T > for all values of of 
x. Due to the smallness of c only small differences occur 
in comparison with the exact values of T, which arc T = 
for x > or which, for the case x < 0, are determined 
by eq.(j§). Hence the set of eqs.@,@ and @ can be 
used for all values of x(t) leading to an unified treatment 
of the cases x < and x > 0. No significant differences 
are found for the results of these two methods. 

In the numerical investigations of the next section the 
results are presented as functions of the temperature from 
T = .01 to T = 1. The temperature variation is per- 
formed both steps-wise and continuously. Usually the 
step size is AT = .01. For the case of continuous vari- 
ations the temperature is linearly and slowly changed. 
Typically 10 5 relaxation times are taken, to change the 
temperature from T = 0.1 to T = 1. 

For runs at a fixed temperature with N = 100 and 
with N — 225, several minutes and a couple of hours are 
needed, respectively. It takes, however, several hours for 
N = 100 and several days for N = 225 when the tem- 
perature is continuously changed in the entire interval. 
At low temperature T < .1 the numerical integration is 
more time consuming than at higher temperatures. 



III. NUMERICAL RESULTS 
A. Finite size effects 

Let us first discuss two trivial finite size effects. The 
largest eigenvalue J max of the interaction matrix J de- 
termines the spin glass temperature T sg = J max /2. For 
finite N this eigenvalue J max differs from the N — > 00 
value II of 2 . The deviations are of the order of 5% and 
mainly arise near T = 1. 

The second simple effect results at low temperatures. 
The binary distributions = i-ZV -1 / 2 imply a splitting 
of the energy levels of AE = 2N~ 1 / 2 . Thus for tempera- 
tures T <C 27V _1 / 2 only the ground state properties enter 
into the calculation of thermodynamic quantities, which 
causes atypical effects for temperature below T = .1. 
These effects are not discussed in the following jT^J . 

For an analysis of the convergence of the present ap- 
proach with increasing N we investigate systems of the 
size N = 16, 25, 36, 49, 64, 81, 100, 121, 196 and 225. For 
each size with N < 121 one hundred different realiza- 
tions of the bonds are investigated. For the systems with 
N = 196 and with N = 225 thirty bond configurations 
are examined in each case. The temperature is fixed at 
T = .5 and all calculations are performed in zero field. 
Starting with random initial values the system relaxes 
always to a fixpoint of the equation of motion that is 
different from the paramagnetic solution m, = 0. 

The fixpoint solutions are classified in two classes. In 
the majority class (containing about 85% of all solutions), 
the value of T is finite and the solutions are located in the 
region x < 0. The solutions of the minority class satisfy 
x > and T = 0. Both classes are already known in 
literature |L2], [l3| and related to the cases where solutions 
of the original TAP equations are not found or are found 

s 

The minority class corresponds to the case were the 
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FIG. 3: Properties of the state of lowest free energy : 

Temperature dependence of \i (full line), of xf^ 1 = — 92) 
(dashed dotted line), of V (dashed line) and of 2 = 1 — /3 2 (1 — 
2g2 + qi) (dotted line) for a zero field cooled system with 
TV = 225 (sample I) which is representative for most systems. 



TAP free energy has a local minimum. As shown in |12 , 
|l3| the obtained x values of these solutions tend to zero 
for N — > 00 and the minimum tends to a marginally 
stable saddle point. Our data are in complete agreement 
with these earlier findings and need therefore no further 
discussion. 

We focus the discussion on the majority class. All so- 
lutions of the this class have a finite values for Y and 
satisfy x < 0. Both results are caused by the finite size 
of the systems, as for N — > 00 all solutions of the mod- 
ified TAP equations with these properties are unstable. 
Thus it is expected that both T and x tend to zero when 
approaching the thermodynamic limit. Fig|l] and Fig|^, 
where the averaged T 2 and the averaged x are plotted 
against N~ 2 ^ 3 , show indeed this expected behavior. The 
data are consistent - but not conclusively - with a TV -2 / 3 
dependence for the asymptotic behavior of both quanti- 
ties. 

From the viewpoint of the modified TAP equations 
both classes are similar. The solutions of the minor- 
ity class and of the majority class describe the approach 
to x — from above and from below, respectively. As 
x = represents the boundary of stability in the thermo- 
dynamic limit, this viewpoint clearly shows the asymp- 
totic marginal stability of the TAP solutions in the glassy 
state. 

Turning to temperature dependent effects we will 
present results for three specific systems. From the ma- 
jority class we select sample I and sample II with sizes 
N = 225 and N = 100, respectively. Sample III with a 
size of N = 100 is chosen from the minority class. 

With the aim of finding the states with lowest free en- 
ergy we slowly cool down the systems to T = .01, starting 
at T = 1 with the paramagnetic solution as initial val- 
ues. For sample I in zero field the numerical results for 
x, T, xi and xf = — 92) are presented in Fig|[ The 
spin glass instability is approached from below for nearly 
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FIG. 4: Properties of the state of lowest free energy: 

Temperature dependence of xi (full line), of xf = (3(1 — 92) 
(dashed dotted line), of Y (dashed line) and of 2 = 1 — /3 2 (1 — 
2?2 +94) (dotted line) for a zero field cooled system with 
TV = 100 (sample III) which is representative for only a few 
systems. 



all temperatures, where x < 0, T > and xi x\ 
hold. The results for sample III are different and are 
plotted in Fig^. In this sample whether the instability 
is approached from above or from below is dependent on 
the temperature. 

In Fig.|| and Fig.|| discontinuities of Xh T and x arise 
which result from saddle node bifurcations and which are 
discussed in some detail in Sec. III.D. These jumps are 
expected to disappear in the thermodynamic limit. In 
this case we have x — ► and r — > and thus Xl(~* 
Xi ) wll l become a smooth curve. Quite remarkably the 
figures show that xi ~* holds for T — > 0, which is a 
desired behavior found in other approaches ^ . 

The procedure of slowly cooling down to get the states 
of lowest free energy can at best be justified by empirical 
arguments. Indeed, further states of the samples under 
investigation were found to have in general, but not al- 
ways a higher free energy (compare Sec. III.C). 

To obtain some more evidence for this procedure, we 
apply our approach to systems with zero fields and with 
N = 10, 11, ... 22 and investigate for each N one hundred 
sets of the bonds. In each case the state of lowest free en- 
ergy at T = is determined by slowly cooling down from 
T = 1 and compared with the exact result which can be 
calculated due to the smallness of the systems. We obtain 
agreement in 95% of all cases and in a further 4% of the 
cases the cooling method found one of the, in general de- 
generate, lowest exited states. A significant dependence 
on N is not found. Provided that an extrapolation to 
large N is possible, these results are quite remarkable, 
even from the viewpoint of optimization problems. 
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FIG. 5: Edwards Anderson order parameter of the 
state of lowest free energy : §2 versus temperature T 
in zero field (full line) and in external fields h for sample I 
with TV = 225 (dotted lines) and for sample II with TV = 100 
(dashed lines). At the temperatures where the lines for finite 
fields end the glassy regime is entered and x becomes nega- 
tive. The dots represent the results of the replica approach 
|g for h = 0. 
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FIG. 6: Total magnetization of the state of lowest 
free energy : m — TV" 1 rm versus temperature in external 
fields h for sample I with TV = 225 (full lines) and for sample II 
with TV = 100 (dashed lines). The dashed dotted line shows 
the SK results ending at the AT temperature and the dots 
mark the temperatures where x changes sign. 



B. Properties of the state of lowest free energy 

In this subsection the results of various physical quan- 
tities will be presented for all temperatures below the 
spin glass temperature in the presence of a homogeneous 
external field h. In some cases the results will exclusively 
be given for sample I and in other cases the results for 
sample II will be added to visualize the finite size effects. 

In Fig.|| the Edwards Anderson order parameter q2 is 
plotted versus T. The finite size effects are rather large in 
the paramagnetic regime, but are reasonable in the spin 
glass regime where q2 is nearly independent of the field 
h. The overall behavior of q2 is in agreement with [ p)[ , 
with the numerical studies of jll| and with the replica 
approach p0|. Note that the agreement with the results 
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FIG. 7: Energy density of the state of lowest free 
energy: u — U/N versus temperature T in external fields h 
for sample I with TV = 225. 
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FIG. 8: Entropy density of the state of lowest free 
energy: s = S/N versus temperature T in external fields h 
for sample I with TV = 225. The dots represent the results of 
the replica approach for h = . 



could be improved by rescaling the temperature axis 
to identical spin glass temperatures. 

In Figj^ the T and h dependence of the total homoge- 
nous magnetization m = N~ 1 '^2 i mi is presented. This 
figure also shows the result of the replica theory [jjj above 
the AT line [fTi] being exact in the thermodynamic limit. 
Again the finite size effects are relatively large. The gen- 
eral behavior, however, can clearly be identified. 

The energy density u = UJN and the entropy density 
s = S/N are plotted in Fig and in Fig.|[ respectively. 
Only the results for sample I are given, as the differences 
to sample II are not resolved on the scale of these figures. 
For the calculation of the plotted quantities eq. ( |l0| ) and 
cq.(|ll]) are used. This use implies again an approxima- 
tion, as strictly these asymptotic equations hold only for 
r = 0. Note finally that the results of FigJ| are again 
in agreement to the results of p0| based on the replica 
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FIG. 9: Meta stable states: Distribution W of the excess 
free energy AF — F — F a at T — .2 in zero field for sample I 
(full line) and for sample II (dashed line). 



approach. 



C. Meta stable states 

Let us estimate the expected number of solutions N s 
of the TAP equations for our system sizes. At T = .2, 
according to Bray and Moore 0], this number is ap- 
proximately given by N s m exp(.05 N). In zero magnetic 
field for each solution a further solution can trivially be 
constructed by changing the sign of all mj. Counting 
each pair as one independent solution the estimates for 
the number of independent solutions are 75 and 3850 for 
N = 100 and N = 225, respectively. Note that the basis 
for this estimate are the original and not the modified 
TAP equations and thus these numbers should be used 
with some care for the present work. 

The investigations for finding meta stable sates are per- 
formed in zero external field at T = .2, starting each run 
with different random initial values of the m,. 

In 540 runs we have found 74 independent solutions of 
the modified TAP equations for sample II. In the first 100 
runs nearly one half of them are obtained. Afterwards 
the rate for finding new solutions quickly decreases to 
about 8% and remains approximately constant for later 
runs. Thus certainly not all existing solutions of this 
sample are observed. The statistical results, presented in 
the following, are nearly independent of the solutions of 
later runs. This implies that the observed solutions are 
representative for all solutions of this sample. A further 
argument for the last statement results from the fact that 
the states with the 8 lowest free energies are found in the 
first 100 runs. This preference for the lower states is also 
found in the rate of recurrence of all runs, which indicates 
a large basin of attraction for these states. 

For sample I with N — 225 we have calculated 868 
different states in 1250 runs. For this sample the rate for 
finding new solutions decreases only to a value of about 
50%. Nevertheless it is assumed that the solutions found 
are again representative for all solutions. In this context 
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FIG. 10: Meta stable states: White-weighted overlap 
distribution at T = .2 in zero field for sample I (full line) 
and for sample II (dashed line). The smooth line shows a 
Gaussian fit of the data of sample I. 



it should be added that about 8% of the meta stable 
states of both samples have a x > and belong therefore 
to the minority class. 

Fig.|| shows the normalized distribution W(AF) of the 
excess free energies AF a = F a — F° of the meta stable 
states, where F a and F° denotes the free energies of the 
meta stable states and of the state of lowest free energy, 
respectively. Note that the range of AF a is identical for 
both samples and seems to be independent of N. More- 
over even the distribution itself seems to be independent 
of N. 

Assuming that the findings can be generalized this is 
an interesting result for the 'multi- valley structure' of the 
free energy. It implies that all the possible excess free 
energies AF a remain finite in the thermodynamic limit. 
Note that this conclusion is in agreement with the result 
of jl0) that the 'multi-valley' structure of the TAP free 
energy occurs on a sub-extensive scale. 

Let us introduce the overlap q aa between the solution 
mf and the solution mf in the usual way 



AT 



(19) 



and define the normalized, white-weighted overlap distri- 
bution 



N S (N S 



T) ? '<« 



q aa ). (20) 



where N s denotes the total number of solutions. (Note 
that P w (q) is not the Parisi overlap distribution, as the 
Boltzmann weights are absent.) 

The distribution P w (q) is symmetrical in q and is plot- 
ted for q > as histogram in Fig. 10 for samples I and for 



sample II. Both distributions are similar and thus nearly 
independent of the system size N. The resulting asymp- 
totic distribution seems to be a Gaussian distribution. 
Such a distribution would result if all the directions of 
the meta stable states in the N dimensional space are 
not correlated. 
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FIG. 11: Meta stable states: Excess free energy AF of 
sample I versus temperature T of the six lowest and four fur- 
ther meta stable states which can be reached by physical pro- 
cesses (compare text). The full and the dotted lines correspond 
to x < and x > 0, respectively and the dashed arrows indi- 
cate the observed jumps when changing T . The insert shows 
the details of a hysteresis near a discontinuity. 



Unfortunately the data for the meta stable states 
are not sufficient to answer the interesting question of 
whether the Parisi overlap distribution, a quantity aver- 
aged over all bonds J5y, applies also to single systems. 
Both samples have just three states with non negligible 
Boltzmann weights and thus definitive conclusions are 
impossible. 

The temperature dependence of the meta stable states 
is presented for (at T = .2) the lowest six and some 
higher states of sample I in Fig.|ll|. The branches of the 
latter states can be reached by physical processes. These 
processes consist of a slow field cooling part which starts 
at T — 1 and which is followed by a sudden switching-off 
of the field at T = .2. For the values h = 1.5, h = 1.0, 
h = .5 and h = .1 the processes end at the points A,B,C 
and D of Fig.[ll], respectively. 

All the meta stable states investigated vary smoothly 
in some temperature region which is limited by an up- 
per limiting temperature. Above this temperature the 
states disappear via saddle node bifurcations. For some 
states an additional lower limiting temperature is found 
(compare insert of Fig.pT|) . The system remains in a def- 
inite meta stable state for all slow changes of the tem- 
perature as long as the limiting temperatures are not 



reached. However approaching a limiting temperature by 
slow temperature changes, the system always bifurcates 
into a branch which has a lower free energy, implying hys- 
teresis effects in some cases. Although such results are 
physically expected this behavior is numerically verified 
for all states plotted in FigJTl]. 

Note that in the whole temperature range most of the 
meta stable states have negligible Boltzmann weights. 
Actually only the three lowest states of sample I will 
give non negligible contributions to Boltzmann averaged 
quantities. Nevertheless the system will remain dynam- 
ically in a meta stable state of high free energy for all 
times provided it remains there initially and provided no 
drastic changes of the temperature or the magnetic field 
are performed. 

According to Fig.|ll| the free energy levels sometimes 
cross with changes of the temperature. Even the state 
of lowest free energy can be involved, which implies that 
the state of lowest free energy and the first exited state 
change their roles at the crossing temperature. Thus in 
the notation of this work, to call the state slowly cooled 
from T = 1 the state of lowest free energy is not strictly 
correct. For brevity we will still use this term keeping, 
however, the limitations in mind. 

D. Saddle node bifurcations 

In this subsection the discontinuities are analyzed 
which arise in the temperature dependence of states of 
lowest free energy and of the meta stable states. 

For this purpose the eigenvalues \^ < \^ < \^ . . . 
and the eigenvectors of the inverse susceptibility jl7| ) 

X'V'"* = A (a) M (a) a = 0,1... N -1. (21) 

are introduced for an arbitrary solution rrii at the tem- 
perature T. Approaching a discontinuity the smallest 
eigenvalue tends to zero and vanishes at a critical 
temperature T c . In Fig.^ an example is presented where 
the temperature dependence of A' ' is given by the nu- 
merical solutions of eq. (pi]) . 

For further discussions the staggered magnetizations 
Am'"' are defined 

AfflW = W T ) - mi(T c )} u\ a) {T c ) a = 0, 1 . . . 

(22) 

The numerical results for Ato'- -' and Am' 1 ' are plotted 
in Fig.|l^ for a discontinuity. Within the numerical errors 
Am(°' ~ (T c - T)' 5 and Am^ - (T c - T) hold where 
the latter result is representative for all a ^ 0. 

These temperature variations can analytically be ex- 
plained. For such an approach eq.(|^) is rewritten as 

hi(T 7 m k ) = — In j ^ ™' - ^ Jijrrij + rriiXi (23) 
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FIG. 12: Saddle node bifurcation: Smallest eigenvalue 
A' ' (dotted line), critical staggered magnetization Am' ' 
(full line) and an uncritical staggered magnetization Am' 1 ' 
(dashed line) versus temperature T near the discontinuity at 
T c = .6707 of Fig.O The dashed dotted line shows the lead- 
ing behavior Am^~ (T c - T) 5 near T c . 




FIG. 13: Free energy landscape: Free energy density 
/ = F/N versus distance v from the critical magnetization in 
direction of u' ' for temperatures above and below the critical 
temperature T c — .591 (compare text). Shown is the heating 
process of sample III for the state of lowest free energy. The 
stable solutions are indicated by dots. The full and the dashed 
lines correspond to the regimes x > and x < 0, respectively. 



and expanded at T c to first order in AT = T c — T and to 
second order in Am; = mj(T) — rrii(T c ) which yields 

AiAT = J2 X^^rn, + ]T ^ Am, Am, (24) 



with 



Ai = d T hi(T,m k )\ T and B ijk = d mk x- 1 (T, mk) I 



(25) 



Multiplying eq.(EJ) from the left by the eigenvectors 
u^ a \T c ) the Am" are calculated to 



Am (0) = ±yJ2AW/BWAT + 0(AT 3/2 ) (26) 
and to 



Am («/o) 



where 



B(°)A(«)(T C ) 



A {a) =Y, u< f ] ( T c) A 



-AT + 0{AT 3/2 ) 



and 



B(<*)=Y,uf\T c )B ijk «f(T„)«W(T e ; 

ijk 



(27) 
(28) 

(29) 



The analytic results ( j2q ) and ( pq ) agree with tempera- 
ture dependence found numerically and plotted in FigJl3. 



Moreover these results exhibit all the typical features of 
a saddle-node bifurcation. The two branches exist only 
for AT > or for AT < depending on the signs of A^ 
and of B(°\ Provided the branches exist one of them is 
stable and the other is unstable p2| . 

The staggered magnetizations Am' 8 ' govern the be- 
haviors of the states near the critical temperatures. Thus 
the Am' ' are the order parameters of the local saddle 
node bifurcations. 

Further numerical evidence for such bifurcations is pre- 
sented in Fig.|l^ where the discontinuity of sample III at 
T = .591 for the heating process is considered. In this 
figure the temperature variation of the free energy land- 
scape is investigated by plotting the quantity 



f(T,v) = N- l F{T, mi -> mi {T c ) + vuf\T c )). (30) 



Below a critical temperature a minimum of the free en- 
ergy and hence a stable solution exists in the region x > 
which disappears by a saddle-node bifurcation above the 
critical temperature. Note that the stationary values of 
v determine approximately the temperature dependence 
of F near T c . 

A further interesting detail is contained in Figjl^. Be- 
low the critical temperature / is not semi-convex every- 
where for x > 0, but above the critical temperature / is 
convex even in some regime x < 0. This demonstrates 
explicitly that the border of stability x — is only ap- 
proximative for finite systems. 
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IV. CONCLUSIONS 

In this work we have demonstrated that the modified 
TAP equations are an adequate tool for exploring the 
characteristics of the pure states of the SK model of fi- 
nite sizes. For all realizations of the random bonds ex- 
plicit solutions exhibiting the spin glass features are ob- 
tained. Moreover, the present approach makes it possible 
to analyze the temperature dependence of all quantities 
of physical interest. In contrast to other approaches the 
method presented can simulate a slow cooling of the sys- 
tem, which is of particular importance to find the features 
of the low lying states. 

Apart from these properties for the state of lowest free 
energy we have found some interesting statistical features 



of the meta stable states which seems to be generic for 
every sample and in the author's opinion, it is a challenge 
to find an analytical approach to these features. 

Due to finite size effects most of the results of this 
work are only qualitative. Thus systems of larger size 
should be investigated. This seems to be at least partially 
possible. 

The present work uses a relaxational dynamics as a 
mathematical tool to find the pure states as fixpoints of 
the equations of motion. Recall that the origin of this 
dynamics is phenomenological and thus true dynamic ef- 
fects cannot be described by the present approach. Thus 
the complete dynamical description should be a subject 
of further research. 
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